Code
library(readxl)
library(EpiEstim)
library(ggplot2)
library(dplyr)
library(janitor)
library(tidyr)
library(knitr)
library(lubridate)
df1 <- read_xlsx("data/data.xlsx", sheet = "data")
df1 <- as.data.frame(df1)Ước tính Rt của Sởi năm 2024
Serial interval
- Dữ liệu hiện có không thể tính được Serial interval theo định nghĩa. Do đó, sẽ dựa vào y văn và chọn Serial interval có trung bình là 14.5 ngày và độ lệch chuẩn là 3.25 ngày (Worden et al. 2020).
Ước tính Rt
Clean data
Code
df1 <- df1 %>% clean_names()
df <- df1 %>% rename(dates = ngay_khoi_phat,
ngaysinh = ngay_sinh,
ngaynv = ngay_nhap_vien_kham,
tiemchung = tinh_trang_tiem_chung,
gioi = gioi_tinh)
sum(is.na(df$dates))[1] 173Code
# Lấy ngày nhập viện thay cho ngày khởi phát của các ca missing
df$dates[is.na(df$dates)] <- df$ngaynv[is.na(df$dates)]
sum(is.na(df$dates))[1] 0- Có 173 ca bị missing Ngày khởi phát nên sử dụng Ngày nhập viện thay thế cho những ca này.
Code
df$cd <- ifelse(df$phan_loai_chan_doan == "Loại trừ sởi", NA, df$phan_loai_chan_doan)
df <- df[,c("dates", "ngaysinh", "tiemchung", "cd")]
# Loại những ca Sởi loại trừ
df <- na.omit(df)
df_convert <- df %>% group_by(dates) %>%
summarise(I = n())
df_complete <- df_convert %>%
complete(dates = seq(min(dates), max(dates), by = "day")) %>%
replace_na(list(I = 0))
df_complete$dates <- as.Date(df_complete$dates)Ước tính hệ số lây nhiễm Rt
- Bộ dữ liệu bao gồm những ca Sởi xác định (Sởi lâm sàng + Sởi xét nghiệm dương tính) từ ngày 04/03/2024 đến 14/09/2024. Sử dụng ngày khởi phát để ước tính Rt .
Code
mod_raw <- estimate_R(
incid = df_complete, # Dữ liệu từ 04/03/2024 - end
method = "parametric_si",
config = make_config(
list(
mean_si = 14.5,
std_si = 3.25
)
)
)Default config will estimate R on weekly sliding windows.
To change this change the t_start and t_end arguments. Warning in estimate_R_func(incid = incid, method = method, si_sample = si_sample, : You're estimating R too early in the epidemic to get the desired
posterior CV.Code
df_rt_raw <- mod_raw$R
df_rt_raw$dates <- mod_raw$dates[df_rt_raw$t_end]
df_rt_raw$q1_rt <- df_rt_raw$`Quantile.0.025(R)`
df_rt_raw$q3_rt <- df_rt_raw$`Quantile.0.975(R)`
df_rt_raw$rt <- df_rt_raw$`Mean(R)`Code
library(patchwork)
p_hist_raw <- ggplot(df_complete, aes(x = dates, y = I)) +
geom_histogram(stat = "identity", binwidth = 1, width = 1, fill = "lightgrey", color = "black") +
labs(x = "Day", y = "Incidence") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) +
scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-03-04"), ymd("2024-09-14")))
p_raw <- ggplot(df_rt_raw, aes(x = dates)) +
geom_ribbon(aes(ymin = q1_rt, ymax = q3_rt), fill = "lightcoral", alpha = 0.3) +
geom_line(aes(y = rt), color = "darkorange") +
geom_hline(yintercept = 1, linetype = "dashed", color = "black") +
labs(x = "Day", y = "Rt", title = "04/03/2024 - 14/09/2024") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) +
scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-03-04"), ymd("2024-09-14"))) +
scale_y_continuous(limits = c(0, 20))
p_hist_raw/p_raw
- Khi ước tính Rt từ 04/03/2024 thì có thể là quá sớm vì các ca bệnh xuất hiện lẻ tẻ từ 04/03/2024 đến khoảng 23/05/2024. Số ca bệnh bắt đầu xuất hiện liên tục từ 23/05/2024 nên sẽ ước tính Rt từ ngày này. Do đó, Rt sẽ được ước tính từ 23/05/2024 đến 10/09/2024.
Code
# Lọc dữ liệu từ 23/05/2024 - end
df_filter <- filter(df_complete, dates >= "2024-05-30")Code
mod_1w <- estimate_R(
incid = df_filter,
method = "parametric_si",
config = make_config(
list(
mean_si = 14.5,
std_si = 3.25
)
)
)Default config will estimate R on weekly sliding windows.
To change this change the t_start and t_end arguments. Code
# Lấy dữ liệu ra để vẽ Rt
df_rt_1w <- mod_1w$R
df_rt_1w$dates <- mod_1w$dates[df_rt_1w$t_end]
df_rt_1w$q1_rt <- df_rt_1w$`Quantile.0.025(R)`
df_rt_1w$q3_rt <- df_rt_1w$`Quantile.0.975(R)`
df_rt_1w$rt <- df_rt_1w$`Mean(R)`Code
library(patchwork)
# Biểu đồ đường cong dịch
p_hist <- ggplot(df_filter, aes(x = dates, y = I)) +
geom_histogram(stat = "identity", binwidth = 1, width = 1, fill = "lightgrey", color = "black") +
labs(x = "Day", y = "Incidence") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) +
scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14")))
# Biểu đồ cho 1 tuần
p_1w <- ggplot(df_rt_1w, aes(x = dates)) +
geom_ribbon(aes(ymin = q1_rt, ymax = q3_rt), fill = "lightcoral", alpha = 0.3) +
geom_line(aes(y = rt), color = "darkorange") +
geom_hline(yintercept = 1, linetype = "dashed", color = "black") +
labs(x = "Day", y = "Rt", title = "30/05/2024 - 14/09/2024") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) +
scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14"))) +
scale_y_continuous(limits = c(0, 20))Code
library(patchwork)
(p_hist_raw + p_hist) / (p_raw + p_1w)
- Trong lệnh
esimate_R, chúng ta có thể tùy chọn Sliding window (SW) khác nhau để thể hiện kết quả ước tính Rt. Sliding window là khoảng thời gian chúng ta dùng để xác định Rt, ví dụ nếu chúng ta chọn SW = 7 thì giá trị Rt sẽ được tính bằng trung bình của 7 ngày. Vậy chúng ta nên chọn *SW bằng bao nhiêu? Chúng ta vẽ thử 4 trường hợp của SW lần lượt là: mỗi ngày, mỗi tuần, mỗi 2 tuần và mỗi tháng để lựa chọn.
Code
# SW = 1
t_start <- seq(2, nrow(df_filter)-1)
t_end <- t_start + 1
mod_day <- estimate_R(
incid = df_filter,
method = "parametric_si",
config = make_config(
list(
mean_si = 14.5,
std_si = 3.25,
t_start = t_start,
t_end = t_end
)
)
)
# Lấy dữ liệu ra để vẽ Rt
df_rt_day <- mod_day$R
df_rt_day$dates <- mod_day$dates[df_rt_day$t_end]
df_rt_day$q1_rt <- df_rt_day$`Quantile.0.025(R)`
df_rt_day$q3_rt <- df_rt_day$`Quantile.0.975(R)`
df_rt_day$rt <- df_rt_day$`Mean(R)`Code
# SW = 14
t_start <- seq(2, nrow(df_filter)-13)
t_end <- t_start + 13
mod_2w <- estimate_R(
incid = df_filter,
method = "parametric_si",
config = make_config(
list(
mean_si = 14.5,
std_si = 3.25,
t_start = t_start,
t_end = t_end
)
)
)
# Lấy dữ liệu ra để vẽ Rt
df_rt_2w <- mod_2w$R
df_rt_2w$dates <- mod_2w$dates[df_rt_2w$t_end]
df_rt_2w$q1_rt <- df_rt_2w$`Quantile.0.025(R)`
df_rt_2w$q3_rt <- df_rt_2w$`Quantile.0.975(R)`
df_rt_2w$rt <- df_rt_2w$`Mean(R)`Code
# SW = 28
t_start <- seq(2, nrow(df_filter)-27)
t_end <- t_start + 27
mod_month <- estimate_R(
incid = df_filter,
method = "parametric_si",
config = make_config(
list(
mean_si = 14.5,
std_si = 3.25,
t_start = t_start,
t_end = t_end
)
)
)
# Lấy dữ liệu ra để vẽ Rt
df_rt_month <- mod_month$R
df_rt_month$dates <- mod_month$dates[df_rt_month$t_end]
df_rt_month$q1_rt <- df_rt_month$`Quantile.0.025(R)`
df_rt_month$q3_rt <- df_rt_month$`Quantile.0.975(R)`
df_rt_month$rt <- df_rt_month$`Mean(R)`Code
library(patchwork)
# Biểu đồ đường cong dịch
p_hist <- ggplot(df_filter, aes(x = dates, y = I)) +
geom_histogram(stat = "identity", binwidth = 1, width = 1, fill = "lightgrey", color = "black") +
labs(x = "Day", y = "Incidence") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) +
scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14")))
# Biểu đồ cho 1 tuần
p_1w <- ggplot(df_rt_1w, aes(x = dates)) +
geom_ribbon(aes(ymin = q1_rt, ymax = q3_rt), fill = "lightcoral", alpha = 0.3) +
geom_line(aes(y = rt), color = "darkorange") +
geom_hline(yintercept = 1, linetype = "dashed", color = "black") +
labs(x = "Day", y = "Rt", title = "1 Week") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) +
scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14"))) +
scale_y_continuous(limits = c(0, 8))
# Biểu đồ cho 1 ngày
p_day <- ggplot(df_rt_day, aes(x = dates)) +
geom_ribbon(aes(ymin = q1_rt, ymax = q3_rt), fill = "lightcoral", alpha = 0.3) +
geom_line(aes(y = rt), color = "darkorange") +
geom_hline(yintercept = 1, linetype = "dashed", color = "black") +
labs(x = "Day", y = "Rt", title = "1 Day") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) +
scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14"))) +
scale_y_continuous(limits = c(0, 10))
# Biểu đồ cho 2 tuần
p_2w <- ggplot(df_rt_2w, aes(x = dates)) +
geom_ribbon(aes(ymin = q1_rt, ymax = q3_rt), fill = "lightcoral", alpha = 0.3) +
geom_line(aes(y = rt), color = "darkorange") +
geom_hline(yintercept = 1, linetype = "dashed", color = "black") +
labs(x = "Day", y = "Rt", title = "2 Week") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) +
scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14"))) +
scale_y_continuous(limits = c(0, 10))
# Biểu đồ cho 1 tháng
p_month <- ggplot(df_rt_month, aes(x = dates)) +
geom_ribbon(aes(ymin = q1_rt, ymax = q3_rt), fill = "lightcoral", alpha = 0.3) +
geom_line(aes(y = rt), color = "darkorange") +
geom_hline(yintercept = 1, linetype = "dashed", color = "black") +
labs(x = "Day", y = "Rt", title = "1 Month") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) +
scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14"))) +
scale_y_continuous(limits = c(0, 5))
(p_hist + p_hist)/(p_day + p_1w)
Code
library(patchwork)
(p_hist + p_hist)/(p_2w + p_month)
Rt cao khi số ca bắt đầu tăng đột ngột và liên tục, càng về sau Rt giảm xuống gần bằng 1 nhưng không duy trì được giá trị Rt < 1 chứng tỏ dịch vẫn đang diễn ra. Với SW = 1 tuần thì Rt ~ 14; Sliding window = 2 tuần thì Rt ~ 19; các khoảng Sliding window khác có giá trị Rt thấp hơn rất nhiều vào ngày ước tính đầu tiên. Theo lý thuyết, Rt sẽ cao khoảng thời gian đầu do dịch mới bùng phát mà 2 tuần trước đó không có ca bệnh xuất hiện liên tục. Do đó, chúng ta chọn giữa SW = 7 hoặc SW = 14 để ước tính.
Tuy nhiên, kết quả cho thấy Rt được ước tính với SW là 2 tuần cho kết quả ổn định hơn, giá trị không lên xuống thất thường nên phù hợp để xem xét và đánh giá diễn tích của dịch. Ngoài ra, khoảng tin cậy của Rt trong thời đầu là quá lớn (khi SW là 1 tuần) nên sẽ không đảm bảo tính tin cậy của Rt.
Ước tính Rt với sliding window là 2 tuần
Code
library(patchwork)
p_hist <- ggplot(df_filter, aes(x = dates, y = I)) +
geom_histogram(stat = "identity", binwidth = 1, width = 1, fill = "lightgrey", color = "black") +
labs(x = "Day", y = "Incidence") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 8)) +
scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14")))
p_2w <- ggplot(df_rt_2w, aes(x = dates)) +
geom_ribbon(aes(ymin = q1_rt, ymax = q3_rt), fill = "lightcoral", alpha = 0.3) +
geom_line(aes(y = rt), color = "darkorange") +
geom_hline(yintercept = 1, linetype = "dashed", color = "black") +
labs(x = "Day", y = "Rt", title = "13/06 - 14/09") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8)) +
scale_x_date(date_labels = "%b %d", date_breaks = "1 week", limits = c(ymd("2024-05-30"), ymd("2024-09-14"))) +
scale_y_continuous(breaks = seq(0,10,by = 2))
p_hist/p_2w
- Dữ liệu để ước tính Rt là từ 30/05/2024, nên Rt sẽ được tính từ ngày 13/06/2024 (Do Serial interval có trung bình là 14 ngày). Rt ~ 7.07 vào 08/06/2024, sau đó giảm dần khi dịch bắt đầu xảy ra (số ca mới xuất hiện liên tục) và cho đến ngày 14/09/2024 (khoảng tuần 36) thì Rt vẫn lớn hơn 1 chứng tỏ dịch vẫn đang diễn ra và số ca bệnh có thể vẫn tiếp tục tăng.
Code
library(knitr)
kable(df_rt_2w[,c("dates", "rt")])| dates | rt |
|---|---|
| 2024-06-13 | 7.0698366 |
| 2024-06-14 | 5.6690111 |
| 2024-06-15 | 4.3508564 |
| 2024-06-16 | 3.2821600 |
| 2024-06-17 | 2.4718480 |
| 2024-06-18 | 2.1445391 |
| 2024-06-19 | 1.8320200 |
| 2024-06-20 | 1.7193623 |
| 2024-06-21 | 1.6256121 |
| 2024-06-22 | 1.5422473 |
| 2024-06-23 | 1.5176277 |
| 2024-06-24 | 1.4954698 |
| 2024-06-25 | 1.3944655 |
| 2024-06-26 | 1.4405630 |
| 2024-06-27 | 1.3418515 |
| 2024-06-28 | 1.1143724 |
| 2024-06-29 | 1.0995479 |
| 2024-06-30 | 1.2085188 |
| 2024-07-01 | 1.4736674 |
| 2024-07-02 | 1.5490903 |
| 2024-07-03 | 1.5791076 |
| 2024-07-04 | 1.5973473 |
| 2024-07-05 | 1.5671345 |
| 2024-07-06 | 1.4948245 |
| 2024-07-07 | 1.4237618 |
| 2024-07-08 | 1.3243248 |
| 2024-07-09 | 1.1656958 |
| 2024-07-10 | 1.0137158 |
| 2024-07-11 | 1.0568080 |
| 2024-07-12 | 1.0603317 |
| 2024-07-13 | 1.1164884 |
| 2024-07-14 | 1.0754648 |
| 2024-07-15 | 0.8953696 |
| 2024-07-16 | 0.9462627 |
| 2024-07-17 | 1.0530632 |
| 2024-07-18 | 1.0882127 |
| 2024-07-19 | 1.1613209 |
| 2024-07-20 | 1.2219629 |
| 2024-07-21 | 1.2947699 |
| 2024-07-22 | 1.4031893 |
| 2024-07-23 | 1.4613128 |
| 2024-07-24 | 1.5730353 |
| 2024-07-25 | 1.6793297 |
| 2024-07-26 | 1.7178643 |
| 2024-07-27 | 1.6554286 |
| 2024-07-28 | 1.6616364 |
| 2024-07-29 | 1.7875178 |
| 2024-07-30 | 1.7556702 |
| 2024-07-31 | 1.6092751 |
| 2024-08-01 | 1.4647907 |
| 2024-08-02 | 1.3967257 |
| 2024-08-03 | 1.3536969 |
| 2024-08-04 | 1.3559618 |
| 2024-08-05 | 1.3996821 |
| 2024-08-06 | 1.5009996 |
| 2024-08-07 | 1.6549297 |
| 2024-08-08 | 1.7120074 |
| 2024-08-09 | 1.8062972 |
| 2024-08-10 | 1.8852197 |
| 2024-08-11 | 1.9690807 |
| 2024-08-12 | 1.8912403 |
| 2024-08-13 | 1.9860504 |
| 2024-08-14 | 1.9307622 |
| 2024-08-15 | 1.9485278 |
| 2024-08-16 | 1.9479556 |
| 2024-08-17 | 1.9839644 |
| 2024-08-18 | 1.9756048 |
| 2024-08-19 | 1.9431634 |
| 2024-08-20 | 1.7682855 |
| 2024-08-21 | 1.6166982 |
| 2024-08-22 | 1.5031903 |
| 2024-08-23 | 1.4052603 |
| 2024-08-24 | 1.3527849 |
| 2024-08-25 | 1.3302986 |
| 2024-08-26 | 1.3423795 |
| 2024-08-27 | 1.2956510 |
| 2024-08-28 | 1.3565252 |
| 2024-08-29 | 1.3408987 |
| 2024-08-30 | 1.2963237 |
| 2024-08-31 | 1.2655504 |
| 2024-09-01 | 1.2486194 |
| 2024-09-02 | 1.2373010 |
| 2024-09-03 | 1.3128349 |
| 2024-09-04 | 1.3294864 |
| 2024-09-05 | 1.3233001 |
| 2024-09-06 | 1.2926711 |
| 2024-09-07 | 1.2373607 |
| 2024-09-08 | 1.1799842 |
| 2024-09-09 | 1.1565229 |
| 2024-09-10 | 1.1126126 |
| 2024-09-11 | 1.0638421 |
| 2024-09-12 | 1.0497415 |
| 2024-09-13 | 1.0363817 |
| 2024-09-14 | 1.0485225 |